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Abstract. The kinetics of the q species pair annihilation reaction {Ai + Aj ~* for 
1 * < i 5: in dimensions is studied by means of analytical considerations and 
Monte Carlo simulations. In the long-time regime the total particle density decays as 
p{t) ~ t^°'. For d = 1 the system segregates into single species domains, yielding a 
different value of a for each q; for a simplified version of the model in one dimension we 
derive a{q) = {q — l)/(2g). Within mean-field theory, applicable in d > 2, segregation 
occurs only for q < 1 + (4/(i). The only physical realisation of this scenario is the two- 
species process {q = 2) in d — 2 and d = 3, governed by an extra local conservation law. 
For d > 2 and q > I + (4/c?) the system remains disordered and its density is shown 
to decay universally with the mean-field power law (a = 1) that also characterises the 
single-species annihilation process A + A ^ 0. 



PACS numbers: 05.40.-a, 82.20.-w 

1. Introduction 

1.1. Motivation 

Thermal equilibrium is a very special situation and nature provides us with 
numerous examples of systems that cannot be described by equilibrium statistical 
thermodynamics. Investigation of simple models is necessary to understand the 
various dynamical phenomena observed in out-of-equilibrium systems, such as collective 
behaviour, phase transitions, or self-ordering [1, 2]. Among these models reaction- 
diffusion processes play a paradigmatic role because of their simple definition on the 
microscopic level and the complex phenomena they may exhibit on the macroscopic 
level [3, 4, 5]. React ion- diffusion processes may describe not only chemical reactions 
[6, 7], but also e.g. phase transitions in ionic conductors [8], epidemic spreading [9], 
or forest fires [10]. The study of such simply defined models is justified by the fact 
that they entail macroscopic phenomena that are universal, i.e. that in the regime of 
asymptotically long times and large distances do not depend on the details of microscopic 
dynamics (for recent reviews, see Refs. [11, 12]). 
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1.2. Annihilation models 

The simplest nontrivial diffusion-limited reaction is the single-species pair annihilation 
process where random walkers (particles of type A) annihilate each other according to 
A+A — > when they meet. Here stands for an inert species that plays no further role 
in the dynamics. The diffusion-limited coagulation process A + A^A displays the same 
asymptotic behaviour [13]. An experimental realisation is the diffusion and annihilation 
of photo-generated excitons on tetramethylammonium manganese trichloride [14]. In 
the long-time regime, for systems with spatial dimension d, the particle density decays 
as p{t) ~ t~^ for d > 2, which is the mean-field result, whereas p(t) ~ t^'^/^ for 
d < 2. Exactly at the upper critical dimension dc — 2 one finds logarithmic corrections 
p{t) ~ t-^\nt. 

These asymptotic power laws were first obtained by means of exact calculations 
for specific models [15]. Their universality is established through a mapping of the 
corresponding classical master equations onto a stochastic field theory [16, 17, 18], 
and its subsequent analysis by means of the dynamical renormalisation group (RG) 
[13, 19]. In the field theory language, the simple annihilation and coagulation processes 
are characterised by the absence of any propagator renormalisation, whence the relevant 
length scale is the diffusion length i{t) ~ t^^"^. The reaction rate enters as a nonlinear 
vertex, whose perturbative renormalisation can be summed to all orders by means of a 
Bethe-Salpeter equation. At long length and time scales, the renormalised reaction rate 
approaches a universal RG fixed point for d < 2. The solution of the Gallan-Symanzik 
equation for the particle density then yields the aforementioned results [19]. 

The two-species annihilation process A -\- B ^ (with no reactions occurring 
between alike particle types) exhibits interesting new phenomena, especially for < 4 
[4, 5, 20]. For a random initial distribution of A and B particles with equal densities 
Pyi(O) = P_b(0), in the course of time segregation into monospecies domains emerges. 
The average domain size is given by the typical diffusion length, and hence grows as 
£{t) ~ t^^^. The reactions happen only in narrow zones between these domains. As 
a consequence the total particle density decays more slowly than in the single-species 
process: Whereas the homogeneous mean-field rate equation solution p{t) ~ t~^ holds 
for d > A, one has p{t) ~ t''^^^ for d < A. The reaction zone width iint{t) also tends to 
infinity according to a power law, imtit) ~ t^'', with, in particular, Xr = 3/8 in d = 1 
[21]. The renormalisation group analysis for this two-species process establishes that the 
upper critical dimension in the RG sense is still dc = 2, i.e., mean-field rate equations 
for the local particle densities augmented with diffusion to account for spatial variations 
provide an apt description for d > 2 [20] . Thus, the borderline dimension for segregation 
^^seg = 4 is located inside the realm of validity of mean-field theory. 

For unequal initial densities of the A and B species [22, 23], the final state will be 
characterised by a nonzero density of the majority species. Above dc — 2 the density 
of the minority species decays exponentially. For d < dc this decay becomes a stretched 
exponential with Inp^i^ ~ —t/lnt foi d — 2 and hip^^^ ~ —f^^"^ for o? < 2. 
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It is worth mentioning that certain special initial conditions with long-range 
correlations arc capable of inhibiting segregation. As an example, consider an alternating 
alignment of particles . . . ABABAB ... in one dimension, and reactions with infinite 
rate. Then any encounter will involve particles of distinct species, leading to their 
annihilation. Yet removing any neighbouring AB pair will conserve the alternating 
alignment. Since ahke particles can therefore never meet, this process is in fact 
equivalent to the single-species reaction A-\- A 0. Segregation is impossible, and 
the total density decays as ~ [24]. 

It is now natural to view the two above models as special cases of a Multispecies 
Annihilation Model (MAM) [25, 26, 27] in which q different species Ai,A2,...,Aq 
perform random walks on a lattice and, when meeting on the same lattice site, may 
react according to 



We will consider that the hopping rates and the reaction rates are species-independent 
{i.e., we assume uniform diffusion constants Di = D and reaction rates Xij — A). For 
equal initial densities Pi(0) = p{0)/q one expects a power law decay 



with a q- and rf-dependcnt exponent that we denote as a{q,d), or, when we are 
interested only in its q dependence, as a{q). In this paper we shall investigate under 
which conditions the MAM may exhibit species segregation, and we will determine the 
value of the exponent a as function of q and d. Thus we fully characterise the long- 
time behaviour of a multicomponent reaction-diffusion system for arbitrary number of 
species, and elucidate the physics behind the different asymptotic scaling laws. Similar 
to other multispecies systems with hard-core constraints [29], the MAM process turns 
out to display very special features in one dimension, owing to the special topological 
limitations of hopping on chains. 

One special limit of the MAM may be analysed immediately. For an infinite number 
of species, q = oo, the probability for two particles of the same type to meet is zero. 
Therefore all encounters lead with rate A to a reaction. Distinguishing different particle 
species then loses its meaning, whence for g = oo we recover again the single-species 
pair process A + A^ in d dimensions [25] . 

We remark that novel behaviour can be expected generically only in the highly 
symmetric case with uniform reaction rates and equal initial densities. In any other 
situation, one would expect that at long times only the least reactive and/or most 
populous species will survive. Once the minority particle species have disappeared, one 
should expect the process to be described by the A + B ^ reaction with unequal 
particle densities [26]. 

1.3. Results 



Ai + Aj^0 , i^j. 




pit) ~ t 



(2) 



We shall present numerical and analytical results concerning the asymptotic decay of the 
total particle density p{t) of the g-species MAM in d dimensions. There are two physical 
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effects that invalidate the homogeneous mean-field rate equation predictions in low 
dimensions, namely reaction rate renormalisation and potentially species segregation. 
Associated with these are two borderline dimensions. First, below the upper critical 
dimension = 2 (independent of q) fluctuations in the dynamics are relevant in the 
RG sense, renormalising the effective reaction rate A. Second, within the mean-field 
approximation, below a g-dependent segregation dimension dgeg given by 

4 

dseg{q) = -^—^ , (3) 

the system splits up into monospecies domains. The behaviour predicted by the 
homogeneous mean-field rate equations (that do not allow for segregation) is observed 
only for d > max(dc, c^seg)- 

The full picture of how the g'-species MAM behaves in d dimensions then emerges 
as follows: 

• For d — 1 the special topological constraints of a hnear chain produce segregation 
with a nonuniversal value of the decay exponent that we determine [26, 27] to be 

a{q,l) = ^, d^l. (4) 

This expression of course reproduces the known one-dimensional exponents 

a(2, 1) = 1/4 and a(oo, 1) = 1/2. 

• For d >2 and g = 2 an extra conservation law in the microscopic dynamics causes 
segregation and 

a{2,d)^^, 2<d<d,eg(2)=4. (5) 

Note that when d = 1 is substituted in Eq. (5), the result coincides with a (2, 1) as 
given by Eq. (4). The local conservation law for the particle density difference in 
two-species annihilation A + B ^ [5] renders the case q — 2 quite special. 

• If one is willing to analytically continue the results to noninteger q, then Eq. (5) is 
a special case of 

a(g,d) = i^^, 2<d<d^g{q), l<q<3. (6) 

For all (g, d) in this range there occurs segregation, although for g 7^ 2 it cannot be 
traced back to a conservation law. 

• For d >2 and q — 3,4, ... ,00 there is no segregation. The density decay therefore 
follows a universal (mean-field) power law, independent of the number of species: 

a{q,d)^l, d>2, q>3. (7) 

In the borderline case d — 2 the power law defined by Eq. (7) is accompanied by a 
logarithmic correction, p{t) ~ Int, precisely as for ^4 -|- A — > 0. 

In one dimension, when segregation occurs, domains of identical particles are 
separated by reaction zones where the annihilation reactions takes place. Let ^int be 
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the typical distance between closest-neighbour pairs A^Aj in a reaction zone. We will 
now assume an infinite reaction rate, so that there is one such pair per zone and i^^i 
is also the zone width. The temporal growth of i-mt may be analysed as follows. The 
density of zones is equal to l/i{t). Hence we may write the rate of decrease of the total 
density as 

dp{t) _ K{t) 

where K{t) is the typical number of annihilations per unit time in a reaction zone. This 
relation, which may be taken as the definition of K,{t), was proposed and exploited by 
Lejrvraz and Redner [21]. The asymptotic time dependence of K{t) is connected by 
Eq. (8) to those of p{t) and £{t). The time needed for the pair AiAj to annihilate is 
~ if^^, so the change of particle density per unit time is dp/dt ~ l/{£if^^). If we now set 
p{t) ~ and £i„t ~ with i ~ t^/^, then Eq. (8) yields the scaling relation [26, 27] 

2q-l 



d=l. (9) 



4? ' 

For the two-species model one recovers ^int ~ t^^^ [21]. 

This paper is organised as follows. In Sec. 2 we recall the known methods and 
results for the q = 2 and q = oo cases in a field-theoretic perspective. We then extend 
these arguments to any value of q. We show why it is difficult in c? < 2 to obtain 
results for general q by renormalisation methods. In Sec. 3 we change our approach 
and consider dimension d = 1 specifically. We introduce a simplified version of the 
MAM, to be referred to as SMAM, in which the stochastic dynamics is replaced by a 
deterministic algorithm. This SMAM is exactly solvable and produces the asymptotic 
density exponent (4). We argue that this result is vahd also for the original stochastic 
MAM. Our analytical findings are supported by our own Monte Carlo simulations, 
presented in Sec. 4, as well as by recent data by Zhong et al. [27]. We show that great 
caution must be exercised when measuring decay exponents, since crossover regimes 
can be long (and depend on the ratio of the reaction rate and the diffusion constant). 
Sec. 5 contains some concluding remarks. Appendix A details the solution of the 
mean-field rate equations. In Appendix B we point out the topological specificities of 
the one-dimensional case by mapping the model onto an anisotropic antiferromagnetic 
Heisenberg spin chain. 



2. Analytical considerations 

2.1. Mean- field theory 

The kinetics of chemical reactions may in first approximation be investigated by 
mean-field theory. The mean-field rate equations neglect any spatial variations of 
the particle densities and thus cannot capture reaction-induced noise and diffusion- 
induced correlations. In low dimensions, however, these effects may govern the long-time 
behaviour of the reaction kinetics, invalidating the mean-field approximation. 
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For the MAM, the mean-field rate equations (with a uniform reaction rate absorbed 
into the time scale) read 

^--jLP^P0^ Z=l,2,...,g. (10) 

For equal initial densities Pi(0) = p*, Eqs. (10) are easily solved to give p{i) — 
rho* /[I + {q — l)p*t]- Above some critical dimension, d > max((ic, dseg) therefore expect 
the MAM decay exponent a to take on the mean-field value Oi{q,d) — 1. 

It is of interest to consider Eqs. (10) when the permutation symmetry of the species 
is broken by the initial conditions, i.e., when the Pi(0) are not all equal. For A; = 1, . . . , p, 
let there be initial values Pfc(O) = p^, having multiplicity n^, where 'Y^k=i''^k — Q- We 
then arrive at a set of p > 1 independent mean-field equations which, after renumbering 
the species and absorbing the rate A into the time scale, may be written as 

^^p- = -Pfc «£Pf - Pfc j , k^l,...,p, (11) 
with the initial conditions now all distinct and ordered such that 

o<pI<pI<...<p;. (12) 

The study of this nonlinear system of equations was initiated by Ben-Avraham and 
Redner [25]. In Appendix A we complete their analysis and determine the full analytic 
solution of the equations (11) and (12) for general p and {n^}. The asymptotic large 
time behaviour may be summarised as follows. 

• For Up > 1 the density Pp of the most abundant species is rip-fold degenerate. It 
decays to zero as t~^. The densities of the other species tend to zero with a faster 
power common to all of them, viz. as 

• For Up = 1 the density of the most abundant species, Pp, is nondegenerate. It tends 
exponentially to a constant PpA* with A* given by 

p-i 

A* = l[{i-p}/p;r. (13) 

e=i 

The densities of all other species tend exponentially to zero for asymptotically large 
times on the common time scale r = 1/ppA*. 



2.2. A Ginzburg-type criterion for the segregation dimension d, 



seg 



The homogeneous mean- field approximation of the preceding subsection 2. 1 is based on 
the hypothesis of there being no particle species segregation. Still within the framework 
of mean-field theory, but allowing now for particle density inhomogeneities, we will 
provide a criterion for the appearance of segregation. To this end, we investigate local 
deviations from the global density averages. When it will appear that these deviations 
are of the same order as the averages, one must conclude that the no-segregation 
hypothesis becomes invalidated. 
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Let there be an infinite system with equal initial densities p* for all q species. We 
will consider a finite subvohimc L'^ of this system. Within such a volume the initial 
densities, due to random fluctuations, will be 



where the Cj are random constants of order unity with zero average. For times t small 
compared to the diffusion time across the volume, t ~ L'^/D, this volume may 
be considered as isolated from the rest of the system. We now assume the mean-field 
approximation to be valid. Generically all q will be different. The PiL{t) then evolve 
in time according to Eq. (11) with all = 1 and initial conditions given by Eq. (14). 
Let Cq be the largest one of the random constants. Hence puit) will tend to zero for all 
i = 1, . . . , g — 1, whereas PqLif) will approach the limit pgL(oo) = p*A*^, with, according 
to Eq. (13), 



Therefore in an isolated volume L'^ only a single species will be left over, with density 
~ The time scale for the approach of this limiting value is tl — 1/ p*A*j^ ~ 

^(9-i)<^/2 This time scale is well within the time that the volume stays isolated at 
the condition that tl ~ ^l, ^-c, if {q — l)d/2 < 2 or equivalently d < dseg, where 
dseg = 4/(g — 1) is the critical segregation dimension announced in Eq. (3). Hence for 
d > dscg diffusion prevents segregation, but for d < dgeg the reactions occur so fast that 
they segregate the system into single-species domains. 

Segregation decelerates the decay of the total density. To determine precisely how 
that happens, we observe that at any time t the total density p{t) is equal to the 
hmit density ^ volume of corresponding size, L{ty, where L{t) ~ {t/Dy^'^. 

Combining this with Eq. (15) we find that asymptotically for t — > oo p{t) ~ t~^i-'i-)d/4:^ 
valid when segregation occurs under mean-field conditions. Hence in this situation we 
obtain the decay exponent a{q, d) = {q — l)d/4 as given in Eq. (6). 

As we will see in subsection 2.3, the mean-field approximation, with or without 
segregation, is valid only for d > dc (and marginally valid in d = d^), where d^ = 2 is 
the g-independent upper critical dimension. It follows that in the mean-field part of the 
{q,d) plane the only physical systems {i.e., with integer q and d) showing segregation 
are the two-species {q — 2) annihilation processes in d — 2 and d — 3 dimensions. We 
cannot exclude, however, that noninteger values of q might someday be given a physical 
interpretation as well. 

It so happens that the q = 2 system satisfies a microscopic conservation law: The 
difference A^^ — Nb between the particle numbers of the two species A and B is strictly 
conserved in each reaction. This law may therefore be viewed as being at the origin 




(14) 



Ai^ll[i-p,L{o)/p,Lm 



i=l 
[q-l 



- HiCq-c) (p*)-(9-l)/2i,-(9-lW2^C'L-(^-lW2. 



(15) 



,i=l 
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of the segregation effect. In any case, for g > 3 and d > 2 we conclude tfiat the 
homogeneous mean-field results of subsection 2.1 apply. 

It is interesting to remark that the derivation of our expression (3) may also be based 
on Ref. [25]. Here, the authors set pi(t) = p(t) + Ci{t)^/]F L^'^^'^ ^ i.e., the time-dependent 
generalization of our Eq. (14), and by expanding the master equation they show that the 
mean-square fluctuations c(t)^ = {{ci{t) — Cj(t))^) obey an ordinary differential equation 
in time [their Eq. (29)]. The solution of that equation, which was also examined very 
recently in Ref. [28], shows that for 1 < ^ < 3 the typical fluctuation c(i)y^L~''/^ 
becomes of the order of the density p{t) on the same time scale tl that appeared above 
in our study of the nonlinear solution; whence we obtain the same result for dgeg- 

2.3. Renormalisation group arguments 

In dimension d < dc — 2 the considerations of the two preceding subsections 2.1 and 
2.2 are no longer valid. Fluctuations in these lower dimensions, where the reactions 
become effectively diffusion-limited, need to be taken into account beyond the Gaussian 
approximation. A systematic treatment of fluctuations in reaction-diffusion problems 
starts from the associated classical master equation, utilises a standard representation in 
terms of bosonic creation and annihilation operators, and therefrom builds a field theory 
action by means of coherent state path integrals, see, e.g., Refs. [16, 17, 13, 19]. This 
action is then further analysed using dynamic renormalisation group (RG) methods, 
and the ensuing scaling laws can be extracted from the Gallan-Symanzik ffow equation. 

For the g-species MAM the ensemble average at time t of any observable F can be 
expressed in terms of functional integrals over 2q fields and as 



total particle density is given by Y^l^i (t>i<t>i)- The statistical weight of any configuration 
is given by the action 



The bilinear part of this action defines q diffusion propagators, whereas the reactions 
introduce q{q — l)/2 two- and four-point annihilation vertices, as depicted in Fig. 1 (a). 
Here, the bare coupling A represents the continuum version of the uniform annihilation 
rate of the discrete model. Straightforward power counting indicates that fluctuations 
become relevant in d < 2 dimensions, for any value of q. Hence we identify dc — 2 SiB 
the upper critical dimension in the RG sense. Thus we need to distinguish three cases: 
(i) For dimensions d > 2 the fluctuation contributions are convergent in the infrared, 
and standard perturbation theory is applicable. Therefore the mean-fleld equations (10), 





(17) 
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J ' '"J J'' ''J J '' J '' 



Figure 1. (a) The propagators and vertices of tlie bosonic field theory. 
(b,c) Renormalisation of the two- and four-point vertices, respectively. 

augmented with diffusive terms, sliould yield the correct asymptotic behaviour, whether 
or not species segregation appears. For values q > 3, the critical segregation dimension is 
dseg < 2, so segregation will not occur. Thus the long-time decay of the particle densities 
will follow the homogeneous mean- field results described in subsection 2.1. Specifically, 
for equal initial densities Pi(0) = p{0)/q the value of the decay exponent is universally 
a = 1, independent of q and d. For the two-species model (g = 2), subsection 2.2 
predicts that segregation happens in dimensions d < 4. As explained in more detail in 
subsection 2.4 below, this is in fact associated with the local conservation of the particle 
density difference pi — p2, which is a very special feature of the two-spccics annihilation 
model. For equal particle densities this leads to the slower decay p{t) ~ t^'^^'^ [4, 5, 20]. 

(ii) At the critical dimension dc — 2 (and below), the full renormalisation group 
machinery is required to further analyse the action (17). As illustrated in Fig. 1(b) 
and (c) , the Feynman diagrams associated with the perturbative renormalisation of the 
annihilation rate assume identical structures for all q. Moreover, the resulting RG flow 
behaves as for the single-species annihilation process A + A ^ 2S (see Refs. [13, 20] for 
the two-species case). In two dimensions, the associated running renormalised coupling 
tends to zero logarithmically for large length- or timescales. Consequently, only the tree- 
level Feynman diagrams involving the initial state enter the computation of the average 
particle densities, cf. subsection 2.4- We may then immediately conclude that for q >3, 
the mean-field results of subsection 2.1 acquire straightforward logarithmic corrections, 
independent of the value of q: If the density of the most abundant particle species is 
nondegenerate {up — 1), it will approach a constant at long times. The approach to 
this asymptotic value, as well as the vanishing of all other species, is characterised by 
a stretched exponential ln[pj(t) — Pi{oo)] ~ — t/lnt. In the degenerate case {up > 1), 
which includes equal initial densities Pi{0) = p{0)/q as a special situation {np = q), 
one obtains p{t) ~ t~^liit. Only the two-species MAM with equal initial densities 
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has distinct features, with p(t) ~ t~^/^, owing to the local conservation law and the 
ensuing particle segregation (sec subsection 2.4)- In Sec. 4, we shall support the above 
predictions by means of Monte Carlo simulations in two and three dimensions. 

(iii) In dimensions d < 2, as explained above, the 'bulk' renormalisation for the 
reaction rate proceeds precisely as for the single-species annihilation model, see Fig. 1(b) 
and (c). However, this does no^ imply that the ensuing asymptotic density decay power 
law holds as well. For now, in order to compute observables such as the averages (pi(t)), 
one must sum all Feynman diagrams that contain contributions from relevant operators 
on the initial time sheet {t — 0). For the g'-species MAM, this turns out to be a 
formidable task. Before we discuss this issue further, let us recall the situation in the 
two-species annihilation model. 

2.4- Effects of the local conservation law for q = 2 

The renormalisation of the annihilation rate A that captures the emergence of particle 
anticorrelations and depletion zones remains present (absent) in dimensions d < dc ~ 2 
{d > 2) and unchanged even if there exists an additional microscopic conservation law. 
Yet physical effects associated with the latter, such as particle segregation, may crucially 
affect the system's leading asymptotic behaviour. Such phenomena, though encoded in 
the initial state, viz. the t — time sheet in the action (17), are in general more difficult 
to analyse in field theory than 'bulk' renormalisations. 

We illustrate this with a brief review of the situation for ^ = 2, following Ref. [20]. 
The two-species pair annihilation process is very special owing to a local conservation 
law (not valid for any other value of q), namely for the particle density difference 
Pi{t) — P2{t) = Pi(0) — P2(0) [5]. As a consequence, the initial state at t = propagates 
into the 'bulk' (t > 0) all the way to t ^ oo, seriously complicating the analytical 
approach. To determine the asymptotic behaviour of the density Lcc and Cardy [20] used 
the following procedure. After shifting the response fields according to 0^ — >■ 0j = 0j — 1, 
new fields ip = (0i + ^2) / V2 and ^ = {(pi — ^2)/ V2 are introduced, along with analogous 
transformations for the response fields tjj and ^. To account for the effect of the initial 
state, specifically for equal initial densities, an effective action is built that contains 
all possible initial terms compatible with the symmetries of the problem. These can 
generically be written as A^"^'"^ ■?/'™(0)^"(0). The most relevant coupling responsible for 
the leading-order decay of the density turns out to be A = A^^'^^ = —A^^''^^ whose 
canonical scaling dimension is 4 — rf. 

For dimensions d > A therefore, the initial state plays no crucial role, and the 
asymptotic behaviour is determined by the homogeneous mean-field rate equations. For 
2 < d < 4, when the 'bulk' parameters in the model do not require any renormalisation, 
a simple effective theory can be constructed. Aside from the new term A^^(O), it 
resembles the mean-field rate equations, supplemented with diffusion terms. Since it is 
a locally conserved mode, the density difference $,(t) propagates diffusively, describing 
the growth of segregated regions in time. This already implies pi{t) ~ t~^^'^ as i — > 00 
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for equal initial densities. The amplitude of this power law is determined by summing 
all tree-level Feynman diagrams that involve the initial parameters p(0) and A, with 
the ultimate result [30, 31] 



In lower dimensions d < 2, however, one needs to apply the full machinery of 
the renormalisation group in order to calculate the (ip) by means of a perturbative 
expansion with respect to the parameter e = 2 — d. Precisely as in the single-species 
pair annihilation process [19], the renormaliscd dimcnsionlcss counterpart of the reaction 
rate A flows towards a stable flxed point of order e. Yet the calculation of observables 
in the framework of this e expansion still necessitates nonperturbative summations over 
all orders of the initial parameters p(0) and A [20]. Whereas this is a straightforward 
task for the single-species reaction, this goal has to date not been achieved beyond 
the leading order in e even for the two-species pair annihilation process. For £ = 0, 
where the renormalised running coupling associated with the annihilation rate vanishes 
logarithmically, this suffices to establish the validity of the result (18) even at the 
upper critical dimension d^ = 2. For d < 2, though, the fleld-theoretic treatment in 
principle leaves three possibilities: Higher-order terms may have no crucial effect at all, 
maintaining the validity of Eq. (18); they may alter the amplitude of the asymptotic 
power law, but leave the exponent a{2,d) = d/i intact; or they might even lead to a 
different asympotic algebraic behaviour. However, this last scenario is not supported 
by simulations in d — 1 [5, 22] and exact results on a specific variant of the two-species 
annihilation model [23]. 

The arguments of Ref. [20] summarised here show that in order to extract the value 
of a{q) in d = 1 within a bosonic field theory, one must perform a sum over all diagrams 
that contain relevant initial state operators. Even for the two-species annihilation model 
this summation appears very difficult, if not impossible. We hence expect the analysis 
to be even more cumbersome for the g-species process. Yet the possibility of species 
segregation must be encoded precisely in these contributions from the initial state. 

To explain the g-dependence of the asymptotic regime, one then needs a field theory 
representation where the special topology of the one-dimensional case becomes manifest. 
An example is provided by the theory of quantum spin chains. In Appendix B, we 
discuss how the annihilation model can be mapped onto an anisotropic antiferromagnetic 
Hciscnbcrg spin chain. The appropriate continuum limit for the low-lying excitation 
yields generalised nonlinear sigma models with topological Chern-Simons terms. For the 
g-species MAM, the different topological sectors for each q will become mixed through 
the annihilation processes. We conjecture that the presence of these topological terms, 
effective only near one dimension, must be responsible for the ^-dependence of the 
density decay exponent in the asymptotic regime. Yet the ensuing formalism is not 
readily amenable to a direct calculation of a(q', 1) either. In view of these difficulties, we 
shall adopt a very different approach to the one-dimensional g-MAM in Sec. 3. In higher 




2 < < 4. 



(18) 
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dimensions d >2 the topological terms do not play any significant role, which supports 
that the RG analysis of section 2.3 predicts the correct MAM asymptotic behaviour. 

2. 5. Role of special initial conditions 

Yet before proceeding, we note that special initial conditions can alter the generic 
asymptotic behaviour of the g-MAM. For instance, if the particles live on a one- 
dimensional lattice, their respective rank can be labeled starting from a given origin. 
Since particles only vanish in pairs, the parity of each particle rank is a locally conserved 
quantity. This immediately implies that if one starts with a correlated configuration such 
as 

. . . Aq^iAgAiA2 . . . Aq^iAqAiA2 . . . , 

the parity of q will dictate the long-time decay of the density. On the one hand, if 
q is odd, two particles of the same species always have an even number of particles 
between each other. Hence these particles can annihilate by pairs and two particles 
of the same species can become nearest-neighbours after some annihilation reactions 
have occurred. Therefore the subtle initial correlations will be wiped out eventually and 
segregation becomes possible, leading to an asymptotic density decay as in the q-MAM. 
(albeit following a potentially large crossover period). On the other hand, if q is even, 
the particles lying between two identical particles cannot entirely annihilate each other. 
In that situation, particles of the same species never meet on adjacent lattice sites, 
species segregation cannot occur, and the density decays as in the single-species pair 
annihilation reaction A + A ^ 0, i.e., with a = 1/2. 

3. Simplified one-dimensional model 

3.1. Simplified Mutual Annihilation Model (SMAM) 

In order to address the one-dimensional case, we now introduce a simplified version of 
the MAM, to be referred to as SMAM, which retains only the bare essentials of the 
original model. We begin by noting that the particle configuration may at any time be 
decomposed into a sequence of domains, each comprising only a single particle species, 
such that adjacent domains contain different species. Owing to the diffusive nature 
of the process the typical domain size increases as l{t) ~ [DtY^'^. as we shall confirm 
numerically in Sec. 4. Combining this with the expectation that p{t) ~ {Dt)~'^, we may 
express the average particle number n{t) = i{t)p{t) in a domain as 

n{t) - (Dt)^-". (19) 

Clearly, unbounded domain growth, hence species segregation, can only occur if a < 1/2. 
The emergence of segregation thus necessarily implies a slower density decay than 
predicted for the homogeneous system {a = 1/2) in one dimension. 

Inspired by Eq. (8) we now introduce the hypothesis that fluctuations in the effective 
annihilation rate K{t), whether in the course of time or between different reaction zones. 
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are not crucial and may be ignored because in one dimension the segregation effect 
will dominate at least the long-time dynamics. This leads us directly to the SMAM 
picture: The particle content of each domain, owing to the annihilation reactions taking 
place at both of its ends, decreases at the constant rate 2K(t). This picture is then 
completed by the rule for what happens when a domain disappears, i.e., gets emptied 
of all its particles: Then, either the left and right neighbouring domains have the same 
species of particles and fuse into a single new domain; or the neighbouring domains have 
different particle species and a new reaction zone appears between them. These two 
possibilities correspond to annihilation and coalescence of domain walls, respectively. 
We already note that it is precisely via these two different scenarios, whose respective 
probabilities depend on the overall number of species, that the decay exponent a{q) 
eventually becomes a function of q. 

The following algorithm makes this explicit. We consider a one-dimensional lattice 
('chain') whose sites, labeled by indices i,j, . . ., represent the domains of the original 
MAM. Below we will indifferently refer to them as 'sites' or as 'domains'. Their 
initial number will be denoted by A^o- To start, we randomly select Nq number 
variables riio, ^20, ■ ■ ■ , JT-iVoO, drawn independently from some probability distribution 
on the positive integers; these variables represent the initial numbers of particles in 
the domains. We then assign initially to each domain a particle type ('colour') in a 
random fashion, except for the restriction that neighbouring domains are not allowed 
the same colour. Yet note that the subsequent algorithm does not require that we assign 
these colours explicitly. The time evolution subsequent to this random initial condition 
proceeds via deterministic iterations. The {k + l)th iteration changes the total number 
of sites from A^^ to A^fc+i and converts the set of values {nik}'^^-^ to {nj^k+i}^^^^ ■ The 
following update consists of four steps: 

(i) The values of all Uik are reduced by 1. 

(ii) All sites i which after step (i) have = are eliminated from the lattice. The 
remaining sites are reconnected in such a way that their ordering along the chain 
remains unchanged. 

(iii) Sites which as a result of step (ii) have become nearest-neighbour sites and which 
have the same colour, are fused into a single site of that colour; the number variables 
of the fusing sites are added together. 

(iv) We write A^fe+i for the remaining number of sites, and relabel these sites by an 
index j that runs from 1 to A^^+i along the chain. 

We began our investigation by actually running this algorithm on the computer, and 
found these simulations to converge rapidly and yield accurate exponent values. These 
advantages were however superseded by our discovery that the above algorithm is 
actually amenable to an asymptotically exact solution. 

In order to keep track of the colours of the sites, in our computer implementations 
of this algorithm we also used a set of colour variables {crik}f=i, where cXik G {1, 2, . . . , g} 
represents the colour of the ith site after the A;th iteration. In our discussion below of the 
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analytic solution of the algorithm our dealing with these variables will be implicit. The 
sole quantity we need to know is the probability, to be denoted rjr, that two domains 

separated by r other domains have the same colour, i.e., that aik = (Tj+r+i,A;- Obviously, 
?7o = and rji = l/{q — 1). Since the two variables a^k and ai^r+i,k can be equal only if 
o'ik 7^ Ci+r,A;, and in that case in fact coincide with probability l/{q — 1), it follows that 
rjr obeys the recursion relation rjr = {1 — '>]r-i)/{(l — 1), whence 

ry, = [l-(l-g)-n/g, r = 0,l,2,.... (20) 

Thus the domain colours exhibit an oscillatory short-ranged correlation. It is easy to 
see that this correlation is in fact invariant under an iteration step. 

After the A;th iteration let = A^^^ denote the average number of particles 

per domain. Let furthermore pk denote the particle density and tk the physical time 
after k iterations. In terms of the data produced by the algorithm, these two quantities 
are given by 

Pk ^ Nkfik tk , . 

po Nofio ' to Ni' ^ ^ 

Here, the second relation holds because the typical domain length i{t) is by hypothesis 
proportional both to the inverse total number of domains, and the square root of time. 
Our purpose will be to determine pk and tk as functions of the iteration index k, and 
subsequently by elimination p{t). 



3.2. SMAM evolution equation 

The key feature of the SMAM responsible for its analytic solvabihty is that at every 
iteration k the number variables nj^, for i = 1,2, . . . , Nk, are uncorrelated. This is 
because they descend from disjoint sets of 'ancestor' variables. As a result this model 
obeys an exact closed set of equations, which we shall now derive and analyse. In 
this subsection we will suppress the iteration index k and will denote, preceding the 
kth iteration, the number of domains containing precisely n particles by the total 
number of domains by A" = Yl'^=i ^n, and the fraction of domains that contain precisely 
n particles by /„ = Mn/N. For the corresponding quantities after the kth iteration we 
shall use primed symbols, i.e., M^, N', and = M'^/N' . By means of appropriate 
combinatorics one may express in terms of the set {fm}m=i before the iteration. We 
now proceed to derive this relationship by investigating the effect of the successive steps 

(i) through (iv) of a single iteration. 

In step (i) there are Mi = Nfi domains that become empty. Let M„ be an auxiliary 
variable indicating the number of domains with n particles (prior to any domain fusion 
processes), and let A^ = Af„ denote the total number of domains left after step 

(ii) of the iteration. Then 

= Mn+i = Nfn+i , n = 1, 2, 3, . . . , 
N ^ N - Ni = N{1- fi) . (22) 
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Next we are interested in the domains that disappear or become created in step 
(iii). Let AA^fusion denote the net change in the total number of domains due to 
this process. This number is easily determined as follows: Since the Nfi eliminated 
domains were randomly placed in the lattice, elementary statistics shows that they 
formed Ns = Nf i{l — fi) sequences of adjacent domains. Their elimination therefore 
led to Ng new nearest-neighbour pairs. The probability Pr for an eliminated sequence to 

contain exactly r domains is = (1 — where r = 1, 2, It follows that two 

newly created nearest-neighbour domains / sites will have equal colours with probability 



p=J2prVr^—:YTh- ^^^^ 

r=l 

We remark in passing that whereas in our earlier paper [26] the possibility of r domains 
disappearing simultaneously was correctly taken into account, the exact expression (23) 
for p was approximated by its low-density, long-time limiting value p — l/{q — 1). The 
exact expression (23) used here actually leads to simpler mathematics, and of course 
does not change the leading power laws; its final consequence is a different amplitude 
in the subleading term of the end result, Eq. (45). 

Since the total number of new nearest-neighbour pairs is N^, the change AA^fugion 
of the total number of domains due to fusion in step (iii) is 

ATVfusion = -N,p = -N . (24) 

Eqs. (22) and (24) allow us to directly relate A^' to N according to 

N'^N + /\N,^,^ = N . (25) 

? - 1 + /i 

However, we have to study step (iii) in greater detail. For we need to know the numbers 
AM„^ioss and AM„^gain of domains of size n that are destroyed and created, respectively, 
during that step. Consider, after step (ii), a specific domain (with arbitrary particle 
number n). This very domain will disappear owing to fusion if its left or its right 
neighbour are of the same colour. The probability p^ that one given neighbour is of 
the same colour is p^ — fi/{q — 1 + /i), namely the product of the probability fi 
of its originally being adjacent to a group of domains of size 1, with the probability 
l/{q — 1 + fi) of this group being adjacent on its other end to a domain of the same 
colour as the specific one we are considering. Since it has two neighbours, the given 
domain under consideration will disappear with probability 2pc — pI- Since there are 
altogether M„ = Nf^+i domains of size n, we find 

AM„,ioss= -M„(2pc-Pe) 

q-1 



1 - 



n=l,2,.... (26) 



Summing this result over n — 1,2,... yields the total number AA^ioss of domains that 
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disappear in step (iii), 



AiVioss = -N{1 - h) 



I + /1 



(27) 



Next we turn to the number AM^ gain of n-particle domains created in step (iii). 
Since such a creation may result from the fusion of s domains, we need to determine the 

corresponding number for each value of s, and subsequently sum over all s = 1, 2, 

After step (ii) there are N — N(l — fi) domains. We may therefore write 

00 

AM„,gain = A^(l - /i) Yl ' ^ = 1> 2, . . . . (28) 

s=2 

Here represents the probability that a specific domain undergoes a fusion with the 
s — 1 domains that follow along the chain, and Rsn is the independent probability that 
this fusion produces a domain of size n. We now determine Xs and Rsn- In order that 
exactly s domains fuse, they must subsequent to step (ii) form a sequence of adjacent 
domains of the same colour, bordered at both ends by domains of a different colour. 
The probability for this arrangement is 

Xs^{l-Po)V-'. (29) 

The domain fractions after step (ii) are /„ = Mn/N = fn+i/0- ~ /i)' "^i^^ ^ = 1> 2, 

The probability that the s fusing domains have sizes mi,m2, . . . ^rus is equal to the 
product fmifm2 ■ ■ ■ frus-i ^^d the probability that this fusion produces a domain of size 
n is in turn equal to this product summed on all sizes of the fusing domains, subject to 
the constraint mi + m2 -\- . . . -\- rus — n. Hence we have 



00 



Rsn — (1 ~ /l) * • • • • • • /ms+1 ^n,mi+...+ms i (30) 

mi=l ms=l 

Using the preceding expressions in Eq. (28) and summing over n — 1,2,... yields 
~ /i(l-/i)(g-l) 

From Eqs. (24), (27), and (31) one may check that AA^fusion = AiVgam + AiVioss, as must 
of course be the case. 

We now consider the recursion 

= M„ + AM„,ioss + AM;,gain , n = 1, 2, . . . , (32) 

and divide both sides of this equation by N' so that the right-hand side becomes equal 
to the fraction f[^. With the aid of several of the preceding equations the right-hand 
side may then be expressed explicitly in terms of the original fractions {/m}m=i- 
thereby finds the SMAM evolution equation 

^^^(i-/i)'(g-i + /i)^"' ^^^^ 
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valid for n — 1,2, . . ., with 

oo oo oo 

'^n = ^ 9l ^ ^ • • • /mi+1 • • • frus+l ^n,mi+...+m„ , (34) 

s=l mi=l ms=l 

and where we have introduced the abbreviation 

Eq. (33), together with (34) and (35), constitutes an exphcit and fully exact recursion 
relation for the evolution of the domain fractions under the deterministic SMAM 
algorithm. One easily checks that it conserves the sum rule Yl'^=ifn = 1- Recall 
that the term with index s in the sum in Eq. (34) represents the creation of a domain of 
n particles by the simultaneous fusion of s domains. The terms with s >3 are obviously 
quite model-specific and one would expect the essential physics to be embodied already 
in the lowest-order nonlinearity, i.e., in the contribution with s — 2. However, in spite 
of the cumbersome appearance of the above recursion, its mathematical analysis turns 
out to be easier when all terms are retained. 



3.3. Asymptotic density decay 

We now restore the iteration indices k and A; -|- 1 to replace the unprimed and primed 
variables, respectively, of the previous subsection. In order to find a solution to 

the SMAM evolution equation (33)- (35) we substitute as an ansatz an exponential 
distribution f^k = efc(l — Cfc)"^^ into the right-hand side of Eq. (33). The recursion 
then reproduces an exponential distribution, fn,k+i = efe-i-i(l — ^^+1)""^, with 

efc+i = ^ ^ — efe . (36) 

q-l + ek 

Hence we have indeed identified a solution. For this solution we may express Eq. (25) 
as 

ta^lH^^^ (37) 

In view of Eq. (21) and the fact that fik — l/cfc, dividing Eq. (37) by (36) gives the 
recursion relation for the particle density pk, 

Pk+i = (1 - ek)pk ■ (38) 

The random initial distribution that we are considering determines the initial condition 
for Eq. (36), viz. eo = (? - l)/q, so /„o = (? - I)/?"- 

Solving Eqs. (36) and (38) explicitly now appears to be quite simple. Prom Eq. (36), 
with the prescribed initial condition, we find 

ek = ^ , (39) 
q + k 

whence, after substitution of this result in Eq. (38), 
Pk _ r(g)r(A: + l) 
po T{k + q) 



(40) 
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Next we may determine tk with the aid of Eq. (21) and the fact that Uk — l/cfc, 



tk ( Po Co 



^0 \Pk^k 



(41) 



_r(g + l) r(A; + l)_ 

The desired time-dependence of the particle density p{t) is then obtained by ehmination 
of the iteration index k from Eqs. (40) and (41) for pk and tk- Wc will be able to carry 
this out explicitly only in an asymptotic expansion for large k. For k ^ q, one easily 
finds from Eqs. (40) and (41) that 

«„ [r(, + i)i= \^ k ^ ) '•■'> 

Setting r = Ctk = [T{q + l)]Hk/to and inverting the expansion (43) gives 

k = tV2« (^1 - i±l r-V2<? + . . . (44) 
Finally, substituting Eq. (44) in (42) and identifying = p{tk) yields 



p{t)^p^V{q) 



2 



(45) 



i.e., the leading-order term of Eq. (45) describes a power-law density decay p{t) ~ 
with the exponent (4); since a{q) < 1/2 we conclude that segregation occurs for all 
2 < g < oo. The correction term in Eq. (45) decays with the power 1/2, which 
corresponds to the single-species decay law in the absence of segregation. Notice that 
indeed the relative amplitude of this term grows with q. In field theory language, 
it represents the effect of reaction rate renormalisation (without segregation) due to 
particle depletion. We note that for large q it will become difficult to numerically 
distinguish this correction from the leading-order term. 

It is not certain that the correction terms that we have identified here for the SMAM 
have the same relevance for the original MAM as we believe the leading-order term 
does. Nevertheless, on physical grounds one would expect the presence of a contribution 
~ Indeed, as we shall see in Sec. 4, Monte Carlo simulation estimates of the MAM 

exponents a{q) in one dimension yield systematically larger values than our prediction 
(4). We interpret this deviation as the influence of an additional term, which 

should also be responsible for the discrepancy between the present theory and earlier 
simulation results [25, 24]. In fact, the detailed analysis of their Monte Carlo data, 
without prior knowledge of our work, led the authors of Ref. [27] to the conclusion 
that the MAM was best described by a superposition of the two power laws present in 
Eq. (45). 



3.4- Reaction zone width 

Let £int be the typical interparticle distance between two unequal particles. In one 
dimension this is also the width of the reaction zone, hypothesised to grow with time 
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as Following the arguments of Ref. [21] as exposed in Sec. 1.2, one is directly 

led to Eq. (9), viz. \r{q) = {2q - l)/(4g). For g = 2 this yields the result A^(2) = 3/8 
derived in Ref. [21] for the two-species system. For q = oo one obtains A(oo) = 1/2, 
which means that the reaction zone width grows as fast as the typical domain size. This 
is possible only if the domains are made up of at most a finite number of particles, in 
accord with the disappearance of segregation as g — > oo and the recovery of q;(oo) = 1/2 
as for the single-species system. 

3.5. Connection to the q-state Potts model 

In the continuum limit colliding domain walls coalesce to a single one with probability 
{q — 2)/{q — 1) and annihilate with the complementary probabihty l/(g — 1). This is 
precisely what happens to domain walls in the dynamics of the one-dimensional zero- 
temperature g-state Potts model. In that context, one may ask for the probability P2{t) 
that after a time t a domain wall that was initially present has not yet disappeared or 
coalesced with another domain wall [32]. In the Potts model one sets P2it) t^^^''\ 
with a g-dcpcndent persistence exponent 6. Although it is known that ^^(2) = 1/2 and 
9{oo) = 3/2, so far an explicit expression for the function 9{q) has not been found (see 
Ref. [33]). 

The same question may be asked for the SMAM. Within this simplified model it 
is easy to find the answer. The calculation is most easily done in the continuum limit. 
Here we state only the result, viz. 9{q) — 1 — 1/q, which is in fact identical to the 
result obtained by Majumdar and Cornell [34] in a mean-field approximation to the 

g-state Potts model. We feel this agreement is not a coincidence, but due to the fact 
that the nonfiuctuating rate hypothesis at the basis of the SMAM actually represents a 
kind of mean-field approximation as well. We conjecture that if appropriate fluctuation 
terms were added to the SMAM evolution equation, then a{q) would remain unchanged, 
whereas 6{q) would acquire nontrivial corrections and deviate from its mean-fleld value. 

4. Monte Ccirlo simulations 

In this section, we report our Monte Carlo simulation results which aim, within the 
numerical accuracy at our disposal, to support the above analytical results concerning 
the exponent a{q,d) in dimensions d = 1, 2, and 3. We begin with a brief description 
of our simulation algorithm. 

First, as appropriate data structure we chose to implement an array-based pseudo- 
binary tree to minimise the time complexity of the following four procedures: selecting 
a random entry, adding new entries, removing entries, and retrieving the number of 
entry matches within the database. In terms of the particle reactions, these processes 
respectively correspond to selecting a random particle, adding new particles at a given 
lattice site, removing particles from a given site, and retrieving the number of particles 
at a given site. This data structure proved to be efficient in single-species reactions. 
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Therefore, it was a natural choice to extend its use to multiple particle species. This 
was accomplished by maintaining separate lists of the particle positions for each species; 
that is, one instance of the tree data structure would store all A particles, one instance 
of the tree data structure would store all B particles, etc. As a consequence of this 
data structure, dropping site exclusion constraints {i.e., allowing only at most one 
particle per site) requires no extra computational effort, and we also initially allow 
particles of different species to coexist on the same site. As the simulation process 
is sufficiently robust to avoid any systematic error in the form of reaction biases, and 
since by diffusion and reaction processes, the probability of multiple occupancy for any 
site diminishes rapidly in time, the asymptotic long-time power laws should not be 
significantly affected by the chosen initial conditions (with the exception of very special 
correlated initial states in one dimension, as will be discussed below). 

Second, the subsequent simulation steps for our implementation of the MAM 
reactions in d dimensions with q reacting species are as follows: 

(i) q (not necessarily empty) sites are selected at random and filled with a particle 
such that exactly one particle of each species is added to the lattice and only one 
particle is added to any of the q sites. 

(ii) Repeat step (i) until the lattice is filled to the desired initial density, usually 1. 

(iii) Increment time by 1/N{t) where N{t) is the total number of particles at this step 
{i.e., apply an asynchronous time update [35, 11]). 

(iv) Pick a random particle, giving equal weight to all particles; more precisely: 

(a) first select the species randomly, weighted by the ratio of the number of particle 
of the given species to the total particle number; 

(b) next, randomly select a particle from this species list, giving equal weight to all 
particles of this species. 

(v) Pick a site randomly from the set of the 2d nearest neighbours of the selected 
particle, observing periodic boundary conditions. 

(vi) If the chosen neighbouring site has no particles of a species different than the 
selected particle, the selected particle simply hops from its original site to the 
neighbouring site. 

(vii) If the chosen neighbouring site contains particles of a species different than the 
selected particle, then one of these particles is randomly selected to be removed 
along with the original particle. 

(viii) Repeat steps (iii) through (vii) until no particles are left or until the overall allotted 
simulation time has elapsed. 

In the following, we will present our simulation results for the g-species MAM with 
g = 2, 3, 4, and 5. The two-species pair annihilation process {q = 2) was investigated 
numerically rather intensively in one [5, 22], two [5, 36] and three dimensions [21]. For 
comparison, we will also show the corresponding data for the single-species annihilation 
process A + A ^ 0. As argued before, this process should describe the g — > oo hmit 
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Figure 2. Local effective decay exponent —a(t) for the one-dimensional MAM 
with g = 2, 3, 4, and 5 (random initial conditions, equal initial particle numbers 
for each species). The data were obtained from 10^ lattice sites, and are averaged 
over 50 runs. For the q = 3 data, we also display the statistical error bars; similar 
accuracies apply to the other graphs as well. The straight lines correspond to the values 
a{q,l) = {q — l)/{2q) as asymptotically predicted by the SMAM. For comparison, 
results for the single-species annihilation process A + A ^ are displayed also, 
corresponding to the limit g — > oo. 



of the g-species MAM in any dimension. Recently, Zhong, Dawkins, and ben-Avraham 
have studied the one-dimensional MAM for g = 3, 4, and 5 by means of a renormahsed 
reaction-cell method for systems that are equivalent to lattices with up to 2^^ (!) sites 
[27]. 

4.I. Segregation in the one- dimensional MAM 

We begin with the MAM in one dimension, where our previous analysis predicts species 
segregation to occur for generic (random) initial placements of the particles along the 
chain. The Monte Carlo simulations were performed on a lattice of 10^ sites with 
periodic boundary conditions, and our data sets were averaged over 50 independent 
runs. For each value of q, the long-time density decay is indeed characterised by different 
exponents a{q,l), see Refs. [26, 27]. In order to avoid any dependence on the details 
of the fitting algorithm and also probe the presence of crossover regimes, we follow the 
standard procedure to measure and plot the 'local' (in time) effective exponent a{t) vs. 
Monte Carlo time t as defined through 

which in the limit t ^ oo should approach a{q, 1). 
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Figure 3. Local effective decay exponent —a{t) for the one-dimensional MAM with 
three particle species {q = 3) with random initial conditions, obtained from averaging 
over 100 runs on 10^ lattice sites. The straight line corresponds to the expected value 
a(3,l) = 1/3. 




Figure 4. One-dimensional MAM for q = 2, 3, and 4. The average size of single- 
species segregated regions grows as t^^'^. 
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Our data as obtained for the one-dimensional MAM with q = 2, 3, 4, and 5 as 
well as for the single- species pair annihilation (corresponding to the limit q oo) are 
depicted in Fig. 2. The horizontal lines indicate the SMAM predictions a(2, 1) = 1/4, 
a{3, 1) = 1/3, a{4, 1) = 3/8, a(5, 1) = 2/5, and a{oo, 1) = 1/2 (from top to bottom). 
While the single-species data agree well with the theoretical prediction, for finite q our 
local effective decay exponent values in the accessible time window are systematically 
lower than the asymptotic SMAM numbers. We interpret these data to reflect the strong 
next-to-leading corrections ~ t~^^'^ to the leading density power law decay, cf. Eq. (45). 
This assertion is in fact borne out by the careful data analysis carried out by Zhong et 
al. for their exceedingly large systems equivalent to lattices of up to 2^® sites. Yet even 
then marked corrections to scaling were still prominent. Notice that the simulation data 
for our system size become rather unreliable for t > 10^. In Fig. 3 we show simulation 
results for the three-species MAM {q = 3) obtained on a larger lattice with 10® sites, 
averaged over 100 runs. The data appear to systematically but slowly approach the 
expected decay exponent value q;(3,1) = 1/3 at long times. 

In addition, we tested a crucial input for our analysis of the SMAM, namely that 
the size of the segregated domains was assumed to grow with time ~ t^^'^. Neglecting 
the effect on the asymptotic regime of the single particle domains, one can look at only 
one realisation, averaging the domain size over all domains in the system to obtain 
almost noiseless data, as shown in Fig. 4 for g = 2, 3, and 4. The measured domain 
growth power law confirms our hypothesis and therefore supports our claim that in the 
long-time regime the SMAM and the MAM become equivalent. 

We also investigated the role of special, highly correlated initial conditions on 
the long-time MAM evolution. As explained in Sec. 2.5, the separation of distinct 
species in an ordered alternating alignment . . . Ag_iAqAiA2 . . . Aq_iAqAiA2 . . . will be 
preserved under the MAM kinetics provided q is even, whence the single-species decay 
with a = 1/2 should be recovered. For odd q, on the other hand, the initial correlations 
disappear after some crossover period, and eventually segregation will occur, with the 
density decay governed by the exponents a{q, 1). These predictions are confirmed by 
our simulation data for q = 2 (on 65536 sites), q = 3 (on 98304 sites), and q = A (on 
65536 sites), each averaged over 20 runs. 

4-2. Simulation results for the MAM in d = 2 and d = 3 dimensions 

For our two-dimensional simulations, we used a 1000 x 1000 sites square lattice with 
periodic boundary conditions in each direction. We have argued above that the single- 
species asymptotic decay law p{t) ~ t^^ \nt should apply for all g > 3 here, in contrast 
with the two-species process for which q;(2,2) = 1/2 as a consequence of the local 
conservation law for the particle density difference pi — p2- In order to test these 
assertions, we plot the local effective density decay exponent a{l) for q — 2, whereas 
for g > 3 we define 

= _dMmiM . (47) 

^ ^ dint ^ ' 
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Figure 5. Local effective decay exponent —o.{t) for the one-dimensional MAM with 
alternating initial particle placement for = 2 (. . . ABABAB ... on 65536 sites), g = 3 
(. . . ABC ABC ... on 98304 sites), and 9 = 4 (. . . ABCDABCD ... on 65536 sites), 

each averaged over 20 nms. For the <? = 3 data, we again display the statistical error 
bars. The straight lines correspond to the theoretical expectations: segregation as for 
random initial conditions only occurs for odd q, while for even q such correlated initial 
states lead to the single-species decay <^ t~^l'^. 



These quantities arc plotted in Fig. 6 for g = 2, 3, 4, 5 and 00 (obtained from the 
single-species reaction), from data averaged over 50 runs. As expected ait) ^ 1 as 
t ^ 00, which confirms the universal density decay for d = 2 and hence the absence of 
species segregation, clearly distinct from the behaviour for q — 2. We consider the data 
reliable up to i = 10^. 

Again, the role of correlated initial states can be investigated. To this end, we ran 
simulations for the two-dimensional two-species annihilation process with checkerboard 
initial particle placement (alternating . . . ABABAB . . . rows) on fairly large square 
lattices with 2048 x 2048 sites, averaging the data over 20 runs. In contrast to its one- 
dimensional counterpart, the stochastic reaction kinetics destroys the initial correlations, 
whence after some crossover period particle segregation emerges and consequently 
pit) ~ at long times as for random initial configurations. 

In three dimensions, the simulations were performed for g = 2, 3, 4, 5, and 00 on 
a 144 X 144 x 144 cubic lattice, and the data averaged over 50 runs. As anticipated 
above, we observe in Fig. 8 the MAM mean-field behaviour with a — 1 for g > 3, 
which demonstrates again the absence of both reaction rate renormalisation and species 
segregation in the system. In contrast, the two-species model should still display particle 
segregation and follow a slower decay with q;(3,2) = 3/4. Notice that while our data 
for q = 2 clearly show behaviour quite distinct from the cases with g > 3, our numerical 
results are markedly off from the expected asymptotic exponent, probably again as a 
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Figure 6. Local effective decay exponent for tlie MAM in d = 2 dimensions for q = 2, 
3, 4, and 5 (random initial conditions, equal initial particle numbers for each species). 
For q>3, the effective exponent a{t) is obtained from p{t)/lnt rather than the density 
itself. The data were obtained on a 1000 x 1000 square lattice, and arc averaged over 
50 runs. For the g = 3 data, we display the statistical error bars. The straight lines 
correspond to the theoretical expectations, p{t) ^ t~^^^ for the two-species case, and 
p{t) ~ t-'^lnt for g > 3. 




Figure 7. Local effective decay exponent —a{t) for the two-species annihilation 
process in d = 2 with checkerboard initial state on 2048 x 2048 sites, averaged over 20 
runs (with statistical error bars). Following an initial crossover period, this correlated 
initial state disappears and segregation ensues with p{t) at long times. 
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Figure 8. Local effective decay exponent —a(t) for the MAM in d = 3 dimensions 
for g = 2, 3, 4, and 5 (random initial conditions, equal initial particle numbers for 
each species). The data were obtained from a 144 x 144 x 144 cubic lattice, and 
are averaged over 50 runs. The straight lines correspond to the values theoretical 
expectations, p{t) ~ t~^/^ for the two-species case, and p{t) ^ for q>3. 



consequence of the competing power law. Beyond t — 10^ the effect of the lattice 
periodicity sets in and our data are no longer useful for extracting the exponent a. 

To summarise, our simulation data confirm the universal g-independent particle 
density decay laws for g > 3 in dimensions d > 2. It is only due to the presence of an 
extra and special conservation law that the two-species model behaves qualitatively 
differently. Correlated initial conditions do not appear to play a crucial role in 
dimensions d>2. 

4-3. Simulation results for two- and three-lane systems 

At last, we were interested in the MAM behaviour on M coupled one-dimensional chains 
(with periodic boundary conditions along and perpendicular to the chains), which we 
refer to as M-lane systems. Within our systematic and statistical errors, and in the 
time window accessible to our simulations, our Monte Carlo data for M = 2 and M = 3 
as depicted in Figs. 9 and 10 yield the same results as for a single one-dimensional 
chain, perhaps with somewhat different initial crossover regimes. Although particles of 
different species may now bypass each other, segregation still occurs for g = 2, 3, 4, and 
5 on such few coupled chains, in contrast to the truly two-dimensional system. 
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Figure 9. Local effective decay exponent —a(t) for the MAM on a two-lane system 
for g = 2, 3, 4, and 5 (random initial conditions, equal initial particle numbers for 
each species). The data were obtained on 2 x 100000 lattice sites, and are averaged 
over 50 runs. For the q = 3 data, we display the statistical error bars. The straight 
lines correspond to the values a{q,l) = {q — l)/{2q) as asymptotically predicted by 
the one-dimensional SMAM. 
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Figure 10. Local effective decay exponent —a{t) for the MAM on a three-lane system 
for g = 2, 3, 4, and 5 (random initial conditions, equal initial particle numbers for each 

species). The data were obtained on 3 x 100000 lattice sites, and are averaged over 
50 runs. For the q = 3 data, we display the statistical error bars. The straight lines 
correspond to the values a{q,l) = {q — l)/{2q) as asymptotically predicted by the 
one-dimensional SMAM. 
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5. Conclusion 

In this paper we have presented the asymptotic long-time behaviour of the particle 
density in the g-species pair annihilation model in d dimensions. Renormalisation group 
arguments show that when q > 3, the total particle density in d > 2 dimensions 
decays according to a universal power law p{t) ~ i"'*, with an exponent a — 1 that 
is independent of q and d. At the upper critical dimension dc — 2, this mean-field result 
is supplemented by logarithmic corrections, viz. p{t) ~ t~^lnt. In dimensions d > 2, 
particle segregation occurs only for the two-species system (g = 2), namely for d < 4, 
resulting in a slower decay with a{2, d) = d/4 [4, 5, 20]. The case q = 2 turns out to be 
quite special owing to the presence of a local conservation law for the difference of the two 
particle densities. Although consistent with the theoretical predictions, our numerical 
investigations in d > 2 could naturally be improved by considering larger systems, 
longer simulation times, and perhaps a more sophisticated algorithm as developed for 
the one-dimensional case in Ref . [27] . 

The topological constraints in one dimension, however, generically induce species 
segregation, and thereby yields a distinct value for the exponent a for each q. A 
simplified, deterministic cellular automaton version of the model allows us to find 
a{q,l) = (q — l)/(2g). This result has been confirmed numerically through Monte 
Carlo simulations, which suggests the equivalence between the original stochastic model 
and the simplified version at least in the accessible time regimes. We have not succeeded 
in calculating the value of a{q,l) using field theory methods, which remains an open 
but difficult task. 

It would also be interesting to study different versions of the MAM in more detail 
[26]. For instance, the g-species cyclic annihilation model (g-CAM) with reactions 
Ai + A2 ^ 0, ^2 + A3 — > 0, . . . , Aq + Ai represents cyclic chemical reactions. 
For low values of q, its properties are easily stated. First of all, the g-CAM and g-MAM 
are straightforwardly equivalent when g = 2 and 3. Then for g = 4, species 1 and 3 
behave in exactly the same manner, and so do species 2 and 4. Hence their respective 
distinction is unnecessary, and the problem is reduced to the 2- MAM, A + B ^ 0. In 
the g = 00 limit, however, the CAM becomes a purely diffusing system with no reaction 
occurring at all. It is then for 5 < g < 00 that the g-CAM and g-MAM can be expected 
to exhibit different behaviour that poses an intriguing problem for further investigation. 
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Appendix A: Solution of the mean-field rate equations (10) and (11) 

We are interested in the solution of the mean- field rate equations (10) with initial 
conditions that break the permutation symmetry between the q particle species. As 
shown in the main text, this leads to a system of p > 1 equations (11) with initial 
conditions Pfe(O) = pj^ ordered according to < p* < P2 < ■ ■ ■ < Pp; where the 
multiplicities Uk obviously satisfy the constraint YTk=i ''^k — Q- (Notice that the case 
p — 2, q — 2, and hence ni = n2 = 1, corresponds to the process A + B 
with unequal initial densities.) We will now obtain the full solution of this system 
of ordinary differential equations. Throughout this appendix, asterisks will refer to 
constants constructed from the initial values. For all 1 < A; 7^ n < p we have p„ 7^ pk 
and may define 

Skn = — ^ , (48) 

Pn - Pk 

and from Eq. (11) it follows that these quantities satisfy 

= -Pn{t) Sknit) . (49) 

at 

Specifically, consider the p — 1 variables Skp {1 < k < p — 1) whose initial values 
Sifcp(O) = sip = pI/ {p* — pI) are ordered according to 

o<«tp<«;p<---<4-i,p- (50) 

Eq. (49) shows that the Skp{t) can only decrease with time but must remain positive. 
Upon dividing two of the equations (49) we find ds^p/ds^p = Sfcp/s^p, which implies that 

^ = |^-4.p, 1<^,^<P-1. (51) 
Sep 

These ratios, of which p — 2 are independent, constitute constants of the motion of the 
nonlinear system of ordinary differential equations (11). 

Before continuing we briefly elaborate on these constants of the motion. Reasoning 
as before one easily shows that in fact all ratios Cken = Skn/ sen (with k, £, n all different) 
are constants of the motion; but they may be expressed in terms of the p — 2 already 
found. There is one more constant of the motion which we do not need explicitly in our 
present calculation but mention for completeness. It is 

c.=pr'flP,-'(fl».»)"""""' (52) 

where n is arbitrary (C^ and C„ with m^n may be expressed through each other with 
the aid of the cun)- Its constancy may be verified by explicit calculation. C„ differs 
from the Cfc£„ in that it is proportional to a nonzero power of the densities. For the 
fully nondegenerate case [p = q and hence all = 1) our Ci reduces to the constant 
of the motion found in Ref . [25] . We emphasise, finally, that when one goes beyond the 
present mean-field approach, these constants of the motion do not correspond to local 
conservation laws, except in the case of q = p = 2, where Ci = l/(pi — P2). 
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As a consequence of Eq. (51), the knowledge of the time evohition of the s^p for a 
single value of k, say Sip(i), suffices. After dividing Eq. (49) for A; = 1 and n — p hy 
(11) for A; = p we obtain 

'^■^ip (53) 

From Eqs. (48) and (51) we have 

Pt-Pp , i^l,...,p-l. (54) 

-Sip + c^^p 

We next use Eq. (54) in (53) to eliminate pi, . . . , pp_i in favour of pp. This results in an 
expression for dsip/dpp in terms of sip and pp. Integration of this expression yields pp 
as a function of Sip. Explicitly, 

Pp = P*p ^i^ip) ' (55) 



where 



'Ip/ l^l \^lp ~ Up , 
Combining (55) with Eq. (49) for A; = 1 and n = p produces the time evolution equation 
for Sip, 

^ = -p;sipF(sip). (57) 

The full solution is now to be obtained as follows. Solving the differential equation 
(57) yields sip(i), which substituted in (55) yields Pp{t). These two results together, 
when used in Eq. (54), produce pe{t) for all £ = 1, ... ,p — 1. For the analysis of (57) it 
is useful to know that the cj^^ are all positive. We remark parenthetically that they are 
in fact ordered according to 1 = Cup > Ci2p > . . . > Ci p_i p > 0. 

We will satisfy ourselves here to extract from Eq. (57) the behaviour of Sip(t) for 
asymptotically large times. In that limit we may set Sip = inside the product, whence 
Eq. (57) reduces to 

^--p;A:,-^sZ, (58) 
where one may verify without great effort that 

= (^tp)-"^'-' 11(1 - pi/p;r = {sip)-'"'^' A* , (59) 

£=l 

with ^4* given by Eq. (13). Now two cases have to be distinguished: 

Case (i): Up — 1. This is the situation where the highest density is nondegenerate. 

Eq. (58) then leads to 

s,p{t)c^C*slpexp{-p;A*t) (60) 

where C* is a constant, and we have factored out Sip. Combining this result with 
Eq. (51) we find that in fact 

s,p(t)^C*s*^exp(-p;A*t) (61) 
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for all /c = 1, ... ,p — 1. Upon substituting the asymptotic solution (60) in Eq. (55) we 
obtain 

p,{t) ~ pp(oo) = p;A* (62) 

Further analysis shows that the approach is exponential on a time scale r = l/p*A*. 
Finally, from the preceding equation and Eqs. (48) and (61), 

Pk{t)c^C*sl^p;A*exp{-p;A*t) 

-C*—^—A*eM-p;A*t), l<k<p-l. (63) 

Pkf Pp 

Conclusion: For Up — 1 the nondegenerate densest species Pp{t) tends exponentially to 
a constant p*A* which depends on all initial densities pi according to Eq. (13). The 
densities of the other species all tend to zero exponentially. The characteristic time scale 
isr=l/(p;^). 

Case (ii): Up > 1. In this situation the highest density is n^-fold degenerate. Eq. (58) 
now leads to 

s,p{t) ^ slp[{np - l)p;A*t]-'/(-^-'^ , (64) 

where we have again factored out a sl^. Combining this with Eq. (51) we find that in 
fact 

sfcp(t) ^ 4p[K - i)p;^*r'/^"^-'^ (65) 

for all /c = 1, ... ,p — 1. Upon substituting the asymptotic solution (64) in Eq. (55) we 
get 

Ppit) ^ p; {s,p{t)/slpr-'A* = [{rip - l)t]-' , (66) 
which is the mean-field behaviour. Finally, from the preceding equation and Eq. (65), 

= ^^[P;^*]-'/^"''-'^ [K - l)t]-i-VK-i) ,l<k<p-l. (67) 
Pp Pk 

Conclusion: For Up > 1 the densest species, which are n^-fold degenerate, tend to zero 
as t~^. The other species tend to zero with a faster power, viz. as 



Appendix B: Mapping onto an antiferromagnetic spin chain in d = 1 

Our goal is to construct a field theory in which the specificity of the one-dimensional case 
becomes manifest. To this end, we recall that the single-species annihilation reaction on 
a lattice can be mapped onto 'quantum' antiferromagnetic XX Z spin chains, provided 
we employ site occupation restrictions, i.e., we allow at most a single particle per lattice 
site X = or 1) [37, 38]. For annihilation processes that asymptotically approach the 
empty, absorbing state, this restriction should not matter for the long-time behaviour. 
Moreover, the RG fixed point for the renormalised reaction rate in ci < 2 dimensions 
actually corresponds to an infinite bare annihilation rate, which precludes multiple site 
occupancy [19]. Instead of using bosonic creation and annihilation operators acting on 
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a Hilbert space that allows arbitrarily many particles per site, we define lowering and 
raising operators a~ and o""*" that operate on the two possible states on each site, 

(7-|0) = , (7+|0) = |1) , (7-|l) = |0) , (7+|l) = , (68) 

and obey the standard anticommutation relations for the spin- 1/2 algebra. 

For the g-MAM, we naturally introduce q operators a^^ which create / destroy a 
particle of species % on site x. It is now a straightforward task to rewrite the classical 

master equation in the form 9t|P(t)) = — /f|P(t)), with a formal state vector \Pif)) = 
E{n.,=o,i} P{{n.,i},t)\{nx,i}) and the pseudo-Hamiltonian H = E.,i(^^:? + E,yi 
with the diffusion part 

Ka = - Y [2nc.,i(n,,+i,i - 1) + (jJ^fT^+i^i + cr-.(7++i_^] , (69) 
and the annihilation contribution 

h7lj = ^ [nx,inx+i,j + nx,jhx+i,i - (J'^i/^x+i,! " '^xj'^x+i^ > (70) 

where we have defined Tix,i — c^jCr^j, and where the summations over i and j extend 
from 1 to q. We remark here that the exclusion constraint is imposed solely on particles 
of the same species. Since the RG fixed point from the 'bulk' renormalisation of A in the 
bosonic field theory describes an effectively infinite microscopic reaction rate, we expect 
that even particles of different species will asymptotically be prevented from occupying 
identical lattice sites at the same (late) time. 

By means of a Jordan- Wigner transformation, one could now proceed to build a 
fermionic field theory from the spin-1/2 pseudo-Hamiltonian H. Indeed, the single- 
species pair annihilation model A + A ^ has been analysed in this framework 
[39, 40, 41, 42]. In order to take into account the topological effects specific to the 
one-dimensional case, we shall follow here the approach of Refs. [43, 44], as described in 

— * 

the Fradkin's textbook [45]. For each spin S, with quantisation axis oriented along Hq 
and eigenstatcs \S, Sz — ±h/2), we introduce a vector n of modulus 1 that defines the 
spin coherent state 

\n) = exp f - ^ ^^^^^ ■ S) \S, Sz = h/2) , with cos^ = no • n. (71) 
V n, smt^ / 

Each vector n has three components (spin coordinates) that we label as n^, n'^, and n^. 
It can be shown that {n\S\n) = nh/2, and that for two different vectors fii and fi2, 



where $(ni, n2, ^o) represents the oriented surface area of the triangle on the unit sphere 
whose vertices are given by the endpoints of the three vectors Hi, ri2, and (fixed) no- 

Through manipulation of these coherent states, the average of an observable can 
be written as in Eq. (16). In discrete time (with the time interval [io,^/] divided into 
N — (tf — to)/St time steps, the hmit 5t — > to be taken in the end) and on the lattice. 
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the ensuing action takes the form 

. L q N ^ 

X=l 1 = 1 k = l ^ 

_ 1 1 + fi.Ahy fixM+i) ^ St {nx,{tkmn.Atk))Y (73) 

Here, the first contribution (with 6t — > 0) is of topological origin, and usually called 
Wess-Zumino term [45]. It is solely responsible for the differences between quantum 
ferromagnets and antiferromagnets, as well as for the remarkable qualitative distinction 
that emerges between half-integer and integer antiferromagnetic spin chains. But notice 
its nonlocal character, since it represents an oriented surface area, yielding disjoint 
sectors each characterised by a certain value of the corresponding topological charge. 
This renders a field theory analysis in terms of a local perturbative expansion futile. 

Let us briefly discuss the single-species pair annihilation model A + A ^ (which 
also describes the g = oo hmit of the q-MAM) in this framework. One obtains exphcitly 

{nx{tk)\H\nx{tk)) = fixitk) ■ nx+i{tk) (74) 



A 

+ 4L 



2<(tfc) - 2<(tfc) <+,(tfc) + i<(t,) <+i(t,) + in:{tk) n'',+,{tk) 



In one dimension, the reactions will be diffusion- limited, i.e., the ratio of the annihilation 
and diffusion rate X/D > 1. Consequently, the isotropic Heisenberg spin coupling in the 
first term is antiferromagnetic (in contrast to a reaction-limited situation, or to a purely 
diffusive model) [39, 40, 41]. Following Ref. [45], we therefore assume that the ground 

— * 

state is represented via the ansatz fix — {—lYrhx + aoix, where ao denotes the lattice 
constant. With this change of variables, we may now study the field fiuctuations around 
the antiferromagnetic ground state. We now take the continuum limit (ao and 
St ^ 0), and integrate out the field i, whereupon after straightforward manipulations 
we arrive at the effective action 5* = -|- Si^^ + S^p, with the harmonic part 



So= dx dt 



— - —mV + (m'^dtm'' - m^dtm") (75) 

ao ao 4ao(A — 2D) 



and the nonlinear contributions 

1 



S,„, = dx dt 



» — » - i6ao(A-2L') 



[m^dtrrf - rrfBtm^f 



X-2D 2 -D 2 
'm'^dtm'^ — m'^dtm^) — - — j- {m'^dtm^ — m^dtm'^) 



4ao(A - 4L>)2 ' ^ ' ' 2ao(A - ADf 

{rrfdtm'^ — m'^dtrrf) {m'^dtm^ — m'^dtm'') 



(76) 



4ao(A-4D)2 

The topological term reads as in Eq. (73), with q = 1. We have furthermore omitted the 
In contribution, since it turns out to be irrelevant. The particle kinetics or, equivalently, 
the spin dynamics can now be examined by diagonalising the quadratic part of the action 
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S'o- To this end, we Fourier-transform Eq. (75) and then determine the eigenvalues 
of the ensuing couphng matrix. Let k and u be the Fourier variables in space and 
time, respectively. We then find that the modes along m'' are static. In the (m^, m'^) 
plane, however, there exist propagating wave solutions with the long-wavelength and 
low-energy dispersion 

Whereas the particle propagation is purely diffusive for A = 0, the presence of 
the annihilation reactions (A > 0) produces an exponential decay of the propagator 
amphtude. Unlike the equivalent bosonic field theory, the full action given by Eqs. (75) 
and (76) is not easily analysed. As in the theory of equilibrium quantum spin chains, 
other methods of analysis such as the Bethe ansatz have proven more fruitful to extract 
the long-time behaviour of such systems. 

We now return to the spin representation of the multi-species annihilation model. 
As evident from Eq. (73), each particle species now carries its own topological charge, 
described by the corresponding Wess-Zumino term. Whereas the diffusion part is still 
diagonal in the species index, 

{n^,iitk)\H\nx,i{tk)y''' = -^n^A'tk) ■ n^+i,i{'tk) , (78) 
the annihilation reactions now couple distinct particle species, 

(n^,i(tfc)|^|n^,i(tfe))^"" = ^Y1 [^^.^(^fc) • ^x+i,j{tk) + <,i(^fc) + Kjitk) 

- 2<*(*^) +^<Atk) <+i,,(4) + ^<,(4) , (79) 

and thus in effect link different topological sectors. This fact further exacerbates the 
difficulties with the spin representation of the g-MAM, and we consider it rather unlikely 
that this formalism will permit a reliable determination of the decay exponent a{q). But 
the presence of the nonlocal topological terms explains why, for integer q at least, the 
q-MAM displays quite different behaviour in d = 1 as compared to higher dimensions. 
For, since we expect the ground state to be characterised by Neel antiferromagnetic 
order, already in two dimensions the contributions of adjacent chains should cancel 
the effects of the Wess-Zumino terms at least at sufficiently long wavelengths [45, 46]. 
For integer q therefore, one would expect special properties, and a dependence of the 
asymptotic scaling regime on the value of q, only in one dimension. 
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